runId = 201;
cnfactType = 6;

%%% Add paths
addpath('../general_funcs');
addpath('../consumerAdoptions');
addpath('../providerAdoptions');
addpath('../providerExits');

%%% This script loads the results of dynamic prediction in the base case and in a set of counterfactual scenarios.
%%% It then produces time series plots to investigate the effect of interventions / what would have happened in alternative scenarios.

%%% Get input / output paths
config = getConfig_jointProcess(runId);

basecaseFolder = sprintf('%s/dynamicPredictions', config.runFolder);

cnfactName = sprintf('cnfact%d', cnfactType);

outputFolder = createFolderIfNotExist(sprintf('%s/cnfacts/%s', config.runFolder, cnfactName));

%%% Load base-case results
load(sprintf('%s/pred_vals.mat', basecaseFolder), 'pred_vals', 'periodLabels');

%%% Load counterfactual results
load(sprintf('%s/cnfacts/%s/cnfact_vals.mat', config.runFolder, cnfactName), 'cnfact_vals');

pred_vals.consumerAdoptions = mean(pred_vals.consumerAdoptions, 3); % T x K
pred_vals.providerAdoptions = mean(pred_vals.providerAdoptions, 3);  % T x K
pred_vals.providerExits     = mean(pred_vals.providerExits, 3);      % T x K
pred_vals.laggedConsumers   = mean(pred_vals.laggedConsumersOrig, 3);   % T x K
pred_vals.laggedProviders   = mean(pred_vals.laggedProviders, 3);    % T x K

cnfact_vals.consumerAdoptions = mean(cnfact_vals.consumerAdoptions, 3); % T x K
cnfact_vals.providerAdoptions = mean(cnfact_vals.providerAdoptions, 3);  % T x K
cnfact_vals.providerExits     = mean(cnfact_vals.providerExits, 3);      % T x K
cnfact_vals.laggedConsumers   = mean(cnfact_vals.laggedConsumersOrig, 3);   % T x K
cnfact_vals.laggedProviders   = mean(cnfact_vals.laggedProviders, 3);    % T x K

pred_vals   = rmfield(pred_vals, 'laggedConsumersOrig');
cnfact_vals = rmfield(cnfact_vals, 'laggedConsumersOrig');


delta.consumerAdoptions = pred_vals.consumerAdoptions - cnfact_vals.consumerAdoptions; % T x K
delta.providerAdoptions = pred_vals.providerAdoptions - cnfact_vals.providerAdoptions; % T x K
delta.providerExits     = pred_vals.providerExits     - cnfact_vals.providerExits; % T x K
delta.laggedConsumers   = pred_vals.laggedConsumers   - cnfact_vals.laggedConsumers; % T x K
delta.laggedProviders   = pred_vals.laggedProviders   - cnfact_vals.laggedProviders; % T x K



%%%%% Compute Effect in IDF vs  Effect outside IDF
% Get thing that associate each kk to a region
indivRunFolders.providerAdoptions  = config.runFolder_providerAdoptions;
indivRunFolders.providerExits      = config.runFolder_providerExits;
indivRunFolders.consumerAdoptions = config.runFolder_consumerAdoptions;
data = load_data_fullmodel(indivRunFolders);

FE_canton = data.consumerAdoptions.Xparts{5}.X_FEs;
is_insideIDF = find(FE_canton == 1);
is_outsideIDF = find(FE_canton ~= 1);


deltaIDF.consumerAdoptions = sum(delta.consumerAdoptions(:,is_insideIDF),2);
deltaIDF.providerAdoptions = sum(delta.providerAdoptions(:,is_insideIDF),2);
deltaIDF.providerExits     = sum(delta.providerExits(:,is_insideIDF),2);
deltaIDF.laggedConsumers   = sum(delta.laggedConsumers(:,is_insideIDF),2);
deltaIDF.laggedProviders   = sum(delta.laggedProviders(:,is_insideIDF),2);

deltaNotIDF.consumerAdoptions = sum(delta.consumerAdoptions(:,is_outsideIDF),2);
deltaNotIDF.providerAdoptions = sum(delta.providerAdoptions(:,is_outsideIDF),2);
deltaNotIDF.providerExits     = sum(delta.providerExits(:,is_outsideIDF),2);
deltaNotIDF.laggedConsumers   = sum(delta.laggedConsumers(:,is_outsideIDF),2);
deltaNotIDF.laggedProviders   = sum(delta.laggedProviders(:,is_outsideIDF),2);

% Put that into a table
M = array2table([deltaIDF.consumerAdoptions, deltaNotIDF.consumerAdoptions, deltaIDF.providerAdoptions, deltaNotIDF.providerAdoptions]);
M.Properties.VariableNames = {'IDF_consumerAdoptionsDelta', 'notIDF_consumerAdoptionsDelta', 'IDF_providerAdoptionsDelta', 'notIDF_providerAdoptionsDelta'};
M2 = cell2table(periodLabels);
M = [M2 M];
writetable(M, sprintf('%s/AdvData.csv', outputFolder));
